* AIR PRODUCT PROJECT
* STEP 1 MODELING: MILP
* LEI ZHANG
* 08/23/2013

SETS
    i   plants      /P1*P2/
    j   customers   /C1*C4/
    k   products    /LOX, LNI, GOX, GNI/
    t   time_points /T1*T4/;

SCALARS
    CU upper bound of expension /100.0/

PARAMETERS

* distance

    PPX(i)  plants_position_x   /P1 100.0, P2 500.0/
    PPY(i)  plants_position_y   /P1 200.0, P2 500.0/
    PCX(j)  customers_position_x    /C1 100.0, C2 300.0, C3 200.0, C4 500.0/
    PCY(j)  customers_position_y    /C1 100.0, C2 100.0, C3 300.0, C4 200.0/
    D(i, j);
    
    D(i, j) = sqrt( (PPX(i) - PCX(j)) * (PPX(i) - PCX(j)) + (PPY(i) - PCY(j)) * (PPY(i) - PCY(j)) );

TABLE DM(j, k, t) customers demand
        T1      T2      T3      T4 
C1.LOX  15.0    16.5    18.2    20.0
C1.LNI  35.0    38.5    42.4    46.6
C1.GOX  15.0    16.5    18.2    20.0
C1.GNI  35.0    38.5    42.4    46.6
C2.LOX  10.0    11.0    12.1    13.3
C2.LNI  30.0    33.0    37.4    41.1
C2.GOX  10.0    11.0    12.1    13.3
C2.GNI  30.0    33.0    37.4    41.1
C3.LOX   5.0     5.5     6.1     6.7
C3.LNI  15.0    16.5    18.2    20.0
C3.GOX   5.0     5.5     6.1     6.7
C3.GNI  15.0    16.5    18.2    20.0
C4.LOX  10.0    11.0    12.1    13.3
C4.LNI  35.0    38.5    42.4    46.6
C4.GOX  10.0    11.0    12.1    13.3
C4.GNI  35.0    38.5    42.4    46.6;

* cost

TABLE ALPHA1(i, k, t) production cost
        T1      T2      T3      T4 
P1.LOX  0.6     0.6     0.6     0.6
P1.LNI  0.6     0.6     0.6     0.6
P1.GOX  0.7     0.7     0.7     0.7
P1.GNI  0.7     0.7     0.7     0.7
P2.LOX  0.5     0.5     0.5     0.5
P2.LNI  0.5     0.5     0.5     0.5
P2.GOX  0.6     0.6     0.6     0.6
P2.GNI  0.6     0.6     0.6     0.6;

TABLE BETA1(i, k, t) investment cost
        T1      T2      T3      T4 
P1.LOX  1.1     1.1     1.1     1.1
P1.LNI  1.1     1.1     1.1     1.1
P1.GOX  1.2     1.2     1.2     1.2
P1.GNI  1.2     1.2     1.2     1.2
P2.LOX  1.0     1.0     1.0     1.0
P2.LNI  1.0     1.0     1.0     1.0
P2.GOX  1.1     1.1     1.1     1.1
P2.GNI  1.1     1.1     1.1     1.1;

TABLE GAMMA1(i, k, t) fixed cost
        T1      T2      T3      T4
P1.LOX  10      10      10      10
P1.LNI  10      10      10      10
P1.GOX  11      11      11      11
P1.GNI  11      11      11      11
P2.LOX  10      10      10      10
P2.LNI  10      10      10      10
P2.GOX  11      11      11      11
P2.GNI  11      11      11      11;

TABLE L(k, t) transition cost
        T1      T2      T3      T4
LOX     0.003   0.003   0.003   0.003
LNI     0.003   0.003   0.003   0.003
GOX     0.004   0.004   0.004   0.004
GNI     0.004   0.004   0.004   0.004;

TABLE PRICE(k, t) price
        T1      T2      T3      T4
LOX     3.0     3.0     3.0     3.0
LNI     3.0     3.0     3.0     3.0
GOX     3.0     3.0     3.0     3.0
GNI     3.0     3.0     3.0     3.0;

TABLE PEN(k, t) competitor price
        T1      T2      T3      T4
LOX     3.1     3.1     3.1     3.1
LNI     3.1     3.1     3.1     3.1
GOX     3.1     3.1     3.1     3.1
GNI     3.1     3.1     3.1     3.1;

VARIABLES
* NPV
    npv;

POSITIVE VARIABLES
* product k produced from plant i to customer j in time t
    x(i, j, k, t)
* production of product k for plant i in time t
    c(i, k, t)
* increased production of product k for plant i in time t
    dc(i, k, t)
    
    z(j, k, t);

BINARY VARIABLES
* if plant i produce product k in time t
    y(i, k, t);

EQUATIONS
    obj, dmd(j, k, t), prd(i, k, t), epd(i, k, t), upb(i, k, t);

obj..   npv =E= SUM( (i, j, k, t), PRICE(k, t) * x(i, j, k, t) ) -  ( SUM( (i, k, t), ALPHA1(i, k, t) * c(i, k, t) + BETA1(i, k, t) * dc(i, k, t) + GAMMA1(i, k, t) * y(i, k, t) ) + SUM( (i, j, k, t), L(k, t) * D(i, j) * x(i, j, k, t) ) + SUM( (j, k, t), PEN(k, t) * z(j, k, t) ) );

dmd(j, k, t).. SUM( i, x(i, j, k, t) ) + SUM(i, z(j, k, t)) =E= DM(j, k, t);

prd(i, k, t)..    c(i, k, t) =G= SUM( j, x(i, j, k, t) );

epd(i, k, t)..   c(i, k, t) =E= c(i, k, t-1) + dc(i, k, t);

upb(i, k, t)..   dc(i, k, t) =L= CU * y(i, k, t);

MODEL APLT /ALL/;

SOLVE APLT USING MIP MAXIMIZING npv;
